#ifndef CSPACECHARGE_H
#define CSPACECHARGE_H
#include <fstream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_halfcomplex.h>

#include <vector>
using namespace std;
const double clight=2.9929e8;
const double pi=3.14159;
class CSpaceCharge
{
    public:
	CSpaceCharge();
	~CSpaceCharge();
	void InputBeam(vector<vector<double> >*,double);
	void InGrid(int,int,double,double);
	void InputStruct(double _dFreq,double _dErest=938.2723);
	void setFreq(double _dFreq);
	void setCellLength(double _dLCell);
	void scheff(double);
	void setPhis(double dPhis=0);
    void getvvd(vector<vector<double> >,string,int iflag1=0,int iflag2=0);
    void formGrid(); 
    void setFlag(int);

    private:
	void initField();
	void EFieldMap();
    void flds(double, double, int, double, double, double, double*, double*);
    void eint(double, double* , double*);
	void plim(double*, double* );
    void Interpret();
    void EFieldMap2();
    void Boundary(int);
	void Potential();
	void FFT2D(vector<vector<double> >*);
	void IFFT2D(vector<vector<double> >*);
	void FFT1D(vector<vector<double> >*);
	void IFFT1D(vector<vector<double> >*);
    public:
    //CBeam beam;
    //int stepNum;  //attention
    private:
    double rHalf,zHalf;
	double dQ,dI;         //粒子电量和正负，现在考虑质子,dQ没用.dI is average current in one rf period.
//        double **compData;
    vector<vector<double> >compData,vvdaa,vvder,vvdez;
    int nr,nz;
    int iflag;
    double dr,dz;
    vector<double> vdzm,vdrm,vdzs,vdrs;
	vector<vector<double> > vvdBeam;            //beam中束流粒子坐标分布
	vector<vector<double> > *cord;            //beam中束流粒子坐标分布
	double xbar,ybar,pbar;
	double dLCell,dPhis,dpLeft,dErest;        //cell同步相位和长度,同步能量
	double dFreq,dWaveLength;
	int iParticleNumTotal,ng;//nip表示计算相邻几个书团之间的空间电荷力，暂时仅考虑书团内的＝1
	double *dEnergyK,dBetaS;
    double dt;
    

    //force==============================================
    vector<vector<double> > vvders,vvdezs;
    vector<int> viemax,vismax;
//    vector<double> vdrssq;
    int nip;
    double c1,pl;
    
    void getPotential(string);
    /*
	double dt,dz,zs,dr,hl,longigridlength,rmax,pl;            //网格尺寸，dt脉冲作用时间
	int nr,nz;              //网格信息
	int iflag;
	vector<double> vdrm,vdrs,vdrssq,vdzm,vdzzs,vders,vdezs; //网格坐标，场强,电荷加权值
	vector<vector<double> > vvdBeam;            //beam中束流粒子坐标分布
	vector<vector<double> > *cord;            //beam中束流粒子坐标分布
	vector<int>  vismax,viemax;
	double xbar,ybar,pbar;
	double *dEnergyK;
	double  dBetaS;
	vector<vector<double> > vvder,vvdez;                       //格点电量的加权值
	vector<vector<double> > vvdaa;
	int iBunchNum;                              //计算时考虑几个书团的作用
	double dLCell,dPhis,dpLeft,dErest,dRadius;        //cell同步相位和长度,同步能量
	double dFreq,dWaveLength;
        double ** dDensityK;
        double ** dDensity;
	int iParticleNumTotal,nip,np,ng;//nip表示计算相邻几个书团之间的空间电荷力，暂时仅考虑书团内的＝1
	//np 是计算过程中遍历用到的变量，ng表示好粒子束。不考虑书束流损失的时候等于iParticleNumTotal
	double c1,c2,c3,c4;
	double xOverR,yOverR,crp,czd;  //Interpret和scheff之间的数据传递
    */
};
#endif
